AD-A174  357 


1/1 


PLASMA  THEORY  AND  SIMULATIONS)  CALIFORNIA  UNIV 
BERKELEY  ELECTRONICS  RESEARCH  LAB  C  K  BIRDSRLL  JUN  85 
NB8814-77-C-S578 

F/6  28/9  NL 


UNCLASSIFIED 


Vvsvw 


AD- A 174  357 


A 


5»jh*s 


l  IJI !  t 


s 


g|  ==3 


B:S  I  o 

:i  sm  c:' 

ill  M 

>>Ni  L-  f~- 
iX'  -  —i 

i^B  :  i 


FIRST  &  SECOND  QUARTER  PROGRESS  REPORT  1985 


ON  PLASMA  THEORY  AND  SIMULATION 


DTIC 

ELECTEJ 

NOV  2  1 1986  . 


January  1  to  June  30,  1985 

DC®  Contract  E£-AS03-76-F00034-DE-A'ITB-76Fr53064 
QNR  Contract  N00014-77-C-0578 
Vanan  Gift 

Hughes  Aircraft  Co.  Gift 


U..  ■. .. 


DISTRIBUTION'S')  AILMENT  A 

Approved  for  pui.hr  trleos*; 
Dis t r. V  u  •  i . >;j  T 


WvV 


••  4 


«•;  *«■ 


ELECTRONICS  RESEARCH  LABORATORY 

College  of  Engineering 

University  of  California,  Bcg'kel^,  ^  ^47^0 


41 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Whan  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT  NUMBER  2.  GOVT  ACCESSION  NO. 

ADAH  ¥357 

3.  RECIPIENT’S  CATALOG  NUMBER 

4.  TITLE  C*n4  SuMItlm) 

C/S&  77T^/r  DA/  Cove*  _ 

S.  TYPE  OF  REPORT  A  PERIOD  COVERED 

Prapcsa,  1/1  -  (WO,  1985 

S.  PERFORMING  ORG.  REPORT  NUMBER 

7.  AUTHOROJ 

Professor  C  K.  Hrcfcall 

S.  CONTRACT  OR  GRANT  NUMBER! A) 

ONR  N00014-77-C0578 

9.  PERFORMING  ORGANIZATION  NAME  ANO  A00RES5 

Electronics  Research  Laboratory 

University  of  California 

Berkeley,  CA  94720 

to.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  4  WORK  UNIT  NUMBERS 

Element  No.  61153N,  Project 

Task  Area  RR01-09-01,  Work 

Unit  No.  NR  012-742 

II.  CONTROLLING  OFFICE  NAME  ANO  AOORESS 

ONR  Physics  Division 

12.  REPORT  DATE 

Department  of  the  Navy,  ONR 

Arlington,  VA  22217 

13.  NUMBER  OF  PAGES 

i4.  .MONITORING  AGENCYN  AME  4  AOClRESS!"  dllUrmtl  trom  Controlling  Ottlco) 

IS.  SECURITY  CLASS,  (of  thia  report) 

Unclassified 

15  a.  DECLASSIFICATION/ OOWNGR  AGING 
SCHEDULE 

l«.  DISTRIBUTION  STATEMENT  (ot  t hi,  Rmport) 

Approved  for  puUic  release;  distribution  unKmifrd 

*7.  DISTRIBUTION  STATEMENT  (ot  the  mOetrmct  entered  In  Block  20,  It  different  from  Report) 

is.  SUPPLEMENTARY  NOTES 

Our  group  uses  theory  and  simulation  as  tools  in  order  to  increase  the  understanding  of  plasma  insta¬ 
bilities,  heating,  transport,  plasma-wall  intnra rtinn«  and  large  potentials  in  plasma*  We  also  work 
an  the  improvement  of  sumdation  both  theoretically  and  practically. 

KEY  WO  AOS  fC«i(lniM  on  r*««rM  aid*  II  no*  earner?  and  Identify  by  block  number) 


Research  in  plwna  thecry  and  simulation,  plasma- wall  interactions,  large  pntmtiah  in  ^annat 

^  ASSTAACT  fCaMku*  mm  rawer  mm  sM  ft  rrmammmaay  mad.  Identity  by  block  number) 


See  reverse  side 


DO»j2T»»W73  KDtTiox  on  nov  «s  is  obsolete  Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  FACE  Dmtm  tnlarnfl 


SSCUWITV  CLASSIFICATION  or  THIS  P AOCfW*«n  O atm  Bntantt) 


20.  ABSTRACT 

r 

i  i  <  ' 

,  i,  --  ■ 

General  Plasma  Theory  and  Simdatian 

A.  Oblique  electron  Bernstein  wave  investigations;  linear  response  lints.  Limits  axe  found  far  a 
spatialy  pea  ode  Vlasov  distribution,  as  a  step  to  understanding  the  bounded  magnetized 
model. 

R  Samlatian  of  the  effect  of  large  anqrlitude  RF -waves  on  the  interchange  instability-supporting 
theory.  Find  report  in  preparation. 

C  One-bean  Alfvin  ion-cyclotron  instabilities  of  multiple  ion  distribution  functions.  Final  report  in 
preparation. 

D.  Linear  made  catqjling  in  sbmlations  of  the  AJfvdn  ion-cyclotron  instability.  Final  report  in 
preparation. 

E  Sarudatians  of  the  weak  warn  beam-plasma  instability :  study  of  the  correlation  functions.  Final 
report  in  preparation. 

Plasma-Wail  Physics,  Theory  and  Sismlatian 

A->  Transport  af  particle  and  energy  fluxes  through  the  plasma  sheath  region,  including  ion  reflec¬ 
tion.  He  rime-tn<fapwiHw*  theory  developed  is  wry  nearly  verified  by  the  average  of  the 
tinvwVpiKlHV  amdationa,  with  exceptions. 

R  Plater  magnetron  discharges.  The  fgperimna  is  now  running  and  initiat  testing  is  being  done, 
electron  density,  I-V  characteristics  M  various  pressures  and  B  fields.  The  simula¬ 
tion  "mrirfling  has  aMrti  volume  creation,  electron-neutral  cdhsians  (tested  to  agree  with 
theory)  and  spstialy  varying  BliekL 

C.  --  Particle  sumdatians  af  the  low-alpha  Pierce  diode.  The  final  report  is  published. 

D, '  Thermionic  emission  I-V  paradox.  The  problem  is  described  (Longris  1/J  observations)  and  con¬ 

jectures  are  offered  on  the  causes. 


Code  Development  and  Software  Distribndan 

A.  A  POLY:  A  hybrid  scheme  far  the  solution  af  the  Vlasov  equation.  One  component  is  modelled  by 
Knearimri  fhad  equations  and  the  cither  components)  by  the  full  Vlasov  equations).  Applica¬ 
tion  is  to  E  Abac.  / 


ft#  Cftft 


dtic 

j^ELECTEgW 

V  NOV  2  1 1986  1 


Copy 

^isPEcmy 


SCCuntTY  CLASSIFICATION  OF  this  P*OEf1Wi«ti  □.-«  Enfrtd) 


TABLE  OF  CONTENTS 


SECTION  L  GENERAL  PLASMA  THEORY  AND  SIMULATION 

A.  Oblique  Electron  Demstcan  Wave  Investigations;  linear 
Response  Limits 

&  Smuiatian  of  the  Effect  af  Large  Amplitude  RF- Waves 
an  the  Interchange  Instability-Supparting  Theory 

C  One-Beam  Alfvdn  Icn-Cyclotron  Instabilities  of 
Multiple  Ion  Distribution  Functions 

D.  Linear  Mode  Coupling  in  Simulations  of  the  AtMa  Ian- 
Cydotran  Instability 

* 

E.  Simulations  of  the  Weak  Warm  Beam-Plasma  Instability: 

Study  of  the  Correlation  Functions 

SECTION  Q:  PLASMA-WALL  PHYSICS,  THEORY  AND  SIMULATION 

A  Transport  Particle  and  Energy  Fluxes  through  the 
Hasma  Sheath  Region,  Including  Ion  Reflection 

R  Planar  Magnetron  Discharges 

C  Particle  Sinmlations  of  the  Low- Alpha  Pierce  Diode 

••• 

D.  Thermionic  Emission  I-V  Paradox 

SECTION  HI:  CODE  DEVELOPMENT  AND  SOFTWARE  DISTRIBUTION 

A  POLY:  A  Hybrid  Scheme  for  the  Solution  of  the  Vlasov 
Equation 

SECTION  IV:  JOURNAL  ARTICLES,  REPORTS,  VISITORS,  TALKS 
DISTRIBUTION  LIST 


Major  support  is  tram  DOR 
^Supported  in  part  by  ONR. 

^Supported  in  pot  by  Vuun. 

Supported  in  pert  by  Hughes  Aircraft  Co. 


‘  f  *«•  *  .  ’  *«'  *,*  \*  *  a  *  »  "  -  «"»*■*  *  *  . 


FIRST  AND  SECOND  QUARTER  PROGRESS  REPORT 

ON 

PLASMA  THEORY  AND  SIMULATION 


January  1  to  Jbae  30, 1985 


Our  research  group  uses  both  theory  and  simulation  as  tools  in  order  to  increase  the  understanding  of 
instabilities,  heating,  transport,  plasma-wall  interactions,  and  large  potentials  in  plasmas.  We  also 
work  an  the  ingnvvement  of  simulation,  both  theoretically  and  practically. 


Our  staff  is 


Professor  C.  K.  Birdsall 

Principal  Investigator 

191M 

Cary  Hall 

(643-6631) 

Dr.  Kim  Theilhaber 

Post-Doctorate;  UCB 

187M 

Cary  Hall 

(642-3477) 

Dr.  Han  Roth  (part-time); 

Research  Physicist,  Space  Science  Lab,  UCB 

187M 

Cary  Hail 

(642-1327) 

Dr.  Bruce  Cohen 

L630 

LLNL 

(422-9823) 

Dr.  A.  Bruce  Langdan 

L472 

LLNL 

(422-5444) 

Dr.  William  Nevins 

Adjunct  Lecturers,  UCB;  Physicists  LLNL 

L630 

LLNL 

(422-7032) 

Mr.  Perry  Gray 

Senior  Engineering  Aide 

119ME 

Cary  Hall 

(642-3528) 

Mr.  William  Lawson 

Mr.  Nids  Otani 

Ms.  Lou  Ann  Schwager 

Research  Assistants 

119MD 

Cory  Hall 

(642-1297) 

119ME 

Cory  Hall 

(642-3528) 

June  30,  1985 

DOE  Contract  DE-ASC3-76-F00034-DE-AT03-76ET53064 
ONR  Contract  N00014-77-C-0578 
Varian  Gift 

Hughes  Aircraft  Co.  Gift 


ELECTRONICS  RESEARCH  LABORATORY 

University  at  California 
Berkeley,  California  94720 


SECTION  I:  PLASMA  THEORY  AND  SIMULATION 


A-  Oblique  Electron  Bernstein  Wave  Investigations 
Linear  Response  Limits 

Win.  S.  Lawson 

In  studying  these  waves,  obtaining  limites  to  linear  response  was  found  to  be  very  important. 
Hence,  a  ampler  model  was  chosen.  The  results  are  being  prepared  as  an  ERL  Memorandum  (now 
ERL  Memo  M86/44,  June  3, 1986)  with  title  and  abstract  as  Mows: 


Limits  of  Linear  Response  of  a  Vlasov  Distribution 

by 

William  S.  Lawson 
Abstract 

The  linear  response  of  a  spatially  periodic  Vlasov  plasma  distribution  function  is  computed  to 
second  order  in  the  electric  field  (this  procedure  is  also  given  justification).  The  results  far  a  specific 
electric  field  are  then  compared  with  the  results  of  computer  simulation  far  different  amplitudes  of 
the  electric  field.  The  onset  of  the  deviations  from  linear  theory  as  the  amplitude  increases  are 
correctly  predicted  by  trapping  theory,  indicating  that  trapping  is  responsible  for  limiting  the  validity 
of  Unear  theory. 


B.  Simulation  of  the  Effect  at  Large  Amplitude  RF-Waves  on 
the  Interchange  Instability-Supporting  Theory 

C.  One-Beam  Alfvfn  Ion-Cyclotron  Instabilities  at  Multiple 
Ion  Distribution  Functions 

D.  Linear  Mode  Coupling  in  Simulations  of  the  Aifv£n 
Ion-Cyclotron  Instability 

Niels  F.  Otani 

Final  reports  are  in  preparation  in  all  of  these  areas,  to  be  issued  as  ERL  memoranda,  to  be 
sent  exit  with  the  QPR’s,  when  ready.  The  tides  may  be  altered  slightly  from  those  above.  One 
abstract  Mows. 


E.  Simulations  of  the  Weak  Warm-Beam-Pi  as  ma  Instability: 

Study  of  the  Correlation  Functions 

Dr.  K.  Theilhaber 

This  is  work  begun  by  Dr.  Theilhaber  in  France  and  completed  here,  to  be 
issued  in  an  ERL  Memorandum,  "Numerical  Simulations  of  Turbulent  Trapping  in 
the  Weak  Beam-Plasma  Instability,"  by  K.  Theilhaber,  ERL  Memo  M86/50,  June  5, 


The  AUVen  Ion- Cyclotron  Instability: 
Simulation  Theory  and  Techniques 


Niels  F.  Otani 

Electronics  Research  Laboratory 
University  of  California 
Berkeley,  CA  94720 


A  particle-ion,  fluid-electron  computer  simulation  code  is  used  in  the  study 
of  the  Alfven  ion-cyclotron  (AIC)  instability,  a  parallel-propagating  electromag¬ 
netic  instability  driven  by  temperature  anisotropy  in  the  ion  velocity  distribution 
function.  A  numerical  odd-even  mode  is  suppressed  by  means  of  a  two-timestep 
averaging  method.  Excellent  energy  conservation  is  obtained  by  using  a  method 
similar  to  the  Boris  particle  mover  to  advance  the  transverse  fields.  Linear  growth 
rates  obtained  from  the  code  differ  substantially  from  those  predicted  by  uniform 
Vlasov  theory,  here  derived  using  a  multifluid  model.  Short  wavelengths  in  partic¬ 
ular  show  substantial  growth  rates  when  damping  is  predicted,  and  additionally 
show  strong  linear  mode  coupling.  Positive  growth  rates  are  even  observed  in 
the  case  of  a  Maxwellian  ion  distribution.  Disagreement  is  also  generally  found 
among  short-wavelength  mode  frequencies.  These  inconsistencies  are  resolved  by 
taking  into  consideration  general  grid  and  discrete- particle  effects  of  the  simula¬ 
tion  model.  A  theoretical  study  reveals  a  real,  physical  process  by  which  an  ion 
distribution  may  collisionlessly  relax  via  short-wavelength  AIC  instabilities  acting 
resonantly  on  small  portions  of  the  distribution  function.  This  process  is  combined 
with  a  linear  mode  coupling  theory  and  other  characteristics  of  the  AIC  instability 
to  explain  all  observed  differences.  Nonlinear  short-wavelength  saturation  levels 
are  also  obtained  and  their  relevance  to  other  field-aligned,  electromagnetic  sim¬ 
ulations  is  discussed. 


To  be  issued  as  ERL  Memo  M85/91,  November  25,  1985 


SECTION  H:  PLASMA- WALL  PHYSICS,  THEORY  AND  SIMULATION 


A.  Transport  of  Particle  and  Energy  Fluxes  through  the 
Plasma-Sheath  Region  Including  Ion  Reflection 

Lou  Am  Scbvager 


I.  Introduction 

Ion  reflection  at  a  collector  plate  affects  the  ability  of  the  sheath  to  thermally  insulate 
the  plate  from  the  plasma.  The  atomic  data  for  ion  reflection  and  derivation  of  electrostatic 
potential  drop  through  the  sheath  and  of  the  kinetic  energy  flux  to  the  plate  is  discussed 
in  a  previous  report  (QPR  III, IV  1984).  Through  this  analysis  and  detailed  comparison 
to  simulation  using  PDWl,  we  have  found  that,  for  a  bounded  system  with  a  constant, 
injection  flux,  increasing  the  reflected  ion  current  generates  a  potential  which  decreases  the 
insulating  effect  of  the  sheath.  For  this  study  we  use  an  ion  to  electron  mass  ratio  of  40 
and  ion  reflection  coefficients  of  R  =  0,  20,  and  40%  with  the  same  system  model  as  before. 

Our  simulation  of  a  bounded  system  with  a  constant  injection  source  models  the  con¬ 
stant  generation  of  ions  and  electrons  by  ionization  in  the  central,  collisional  section  of  a 
plasma  device.  Consequently  boundary  effects,  such  as  ion  reflection,  affect  not  only  the 
plasma  characteristics  but  also  the  actual,  equilibrated  density  in  the  system.  Thus  for  a 
system  with  constant,  injected  flux,  the  increase  in  ion  reflection  causes  a  buildup  of  plasma 
density  which  overall  increases  the  kinetic  energy  flux  to  the  plate. 

II.  Theory  Review 

In  the  previous  report  (QPR  HI, IV  1984)  we  describe  the  kinetic  energy  flux  to  the 
plate  as  the  sum  for  all  three  plasma  components  of  the  energy  carried  per  particle  times  the 
particle  flux  at  the  electrically  insulated  boundary.  First  let  us  consider  the  contribution 
of  the  electrons.  In  one  dimension,  on  average  each  incident  electron,  with  a  mean  speed  of 
{v)  =  (2kT/wm)1/2,  carries  an  energy  of  IkT  to  the  plate  by  thermal  effusion.1  Kinetically 


the  mean  speed  of  these  Maxwellian  electrons  is  defined  by 


{v)  =  (2tt vf)  1/2  J v  exp  (1) 

o 

The  density  of  electrons  hitting  the  plate  is  one-half  the  density  determined  by  the  electro¬ 
static  potential  through  the  Boltzmann  equation: 


1  _ _ /  e4>o\ 

ne°  ~  2n°°  eXP  \  kT  ) 


(2) 


A  discussion  of  why  the  electrons  carry  an  energy  of  1  kT  (in  one-dimension)  is  presented 
in  the  appendix  of  the  aforementioned  report.  Hence  electrons  with  an  energy 


Ee0  =  kT 


(3) 


contribute  a  kinetic  energy  flux  of 


QeO  =  FeokT, 


(4) 


where  the  particle  flux  at  the  plate  is 

F*o  =  exp  (”)■  (5) 

Using  the  following  definition,  we  can  derive  this  relation  for  Qe o  m  one-dimension  with 

OO 

QeO  =  noo(2™*)~l/2  J  v3exp(^^^jdv.  (6) 

o 

We  next  evaluate  the  contribution  of  the  primary  and  reflected  ions  to  the  kinetic 
energy  flux  at  the  plate.  Because  we  assume  the  plate  is  electrically  insulated  or  floating, 
the  electron  flux  Fe o  equals  the  primary  ion  flux  F,o  into  the  plate  minus  the  reflected  ion 
flux  Fro  away  from  the  plate.  Using  the  definition  of  ion  reflection  coefficient, 


along  with  the  open  circuit  condition,  we  find  that 


Fi0  =  FM/(l  -  R) 


(8) 


and 

-  -RFeo/il  -  R).  (9) 

The  primary  ions  cany  their  initial  energy  Eioo  on  entering  the  sheath  region  plus  that 
gained  from  the  drop  in  potential  of  —e<t>Q.  Thus  each  primary  ion  imparts  to  the  plate  a 
total  energy  £,a  where 

£i0  =  Zoo  ~  efo.  (10) 


For  simplicity  we  assume  the  plate  absorbs  no  energy  during  reflection.  *  Then  the  reflected 
ions  remove  the  same  energy  as  the  primary  ions  add.  As  shown  in  the  previous  report 
(QPR  m,IV  1984:  Eq.  (18)),  £ioo  >  kT/2  for  cold  ions.  With  the  injection  of  warm  ions, 
T,  =  T,  these  enter  the  sheath  region  with  initial  energy  £««>  equal  to  kT.  With  this 
assumption  then  the  kinetic  energy  flux  for  primary  ions  is 


and  for  reflected  ions  is 


.  -RFepkT  (y  e<^\ 

Qr0~  1  -R  V  kT)' 


In  an  equilibrated  system  with  a  constant,  injected  flux  and  temperature,  the  electron 
flux  incident  to  the  plate  must  be  equal  to  the  injected  flux  Ftoo.  With  the  same  F -o  for 
*  Typically  for  a  low  temperature  plasma,  the  wall  absorbs  about  20%  of  the  incident 
particle  energy.  The  specific  fraction  depends  on  the  mass  ratio  of  the  incident  particle  and 


each  variation  in  reflection  coefficient  R,  we  see  that  Qe 0  is  independent  of  R  as  is  n^. 
Thus  for  the  plate  electron  density  of  Eq.  (2)  to  be  independent  of  R  (and  4>o)  then 
must  be  increasing.  For  the  ions,  net  kinetic  energy  flux,  Q,o  +  Qro,  depends  directly  on 
the  floating  plate  potential.  As  we  increase  R,  this  drives  up  the  magnitude  of  the  potential 
4>o  according  to  the  previously  derived  relation  (QPR  m, IV  1984:  Eq.  (21)),  which  is 


Thus  the  net  ion  contribution  rises  slightly  only  because  of  the  added  energy  that  the  ions 
gain  in  accelerating  through  a  larger  potential  drop. 

Instead  of  fixing  the  injected  flux,  suppose  we  hold  constant  by  appropriately 
lowering  the  injected  fluxes  as  we  increase  R.  Then  the  larger  sheath  drop  will  better 
insulate  the  plate  from  the  hot  electrons,  i.e.  Qe 0  falls  linearly  with  Feoo.  With  fixed, 
ion  kinetic  energy  flux  falls  because  the  exponential  drop  in  plate  density  overshadows  the 
small  gain  in  incident  ion  energy. 

The  choice  of  whether  to  fix  or  Feoo  dramatically  affects  the  dependence  of  Q0 
on  R.  We  show  this  difference  in  the  two  graphs  of  Fig.  1.  The  parameters  used  in  the 
simulation  and  substituted  into  Eqs.(4)  and  (11)-(13)  for  Fig.l  are  F^q  =  500,  M/m  =  40, 
and  e  =  kT  =  0.05.  Potential  is  calculated  with  Eq.  (13).  We  believe  that  holding  the 
injected  flux  constant,  rather  than  n^,  better  models  a  bounded  plasma  device.  From  this 
point  forward  all  statements  of  trends  will  use  the  assumption  that  we  fix  the  injected  flux 
while  varying  R. 

III.  Simulation  Transport  Diagnostics 

Next  we  present  how  PDWl  evaluates  the  above  plasma  characteristics.  Diagnostics 
for  moments  of  the  distribution  function  are  evaluated  as  particles  leave  the  bounded  sys¬ 
tem  within  each  time  step.  Usually  at  each  time  step  we  can  extract  the  velocity  of  all 


> 


the  particles  weighted  to  a  grid  point.  Distributing  the  particles  according  to  velocity  into 
many  velocity  bins  of  width  An  determines  at  that  position  and  time  the  particle  distri¬ 
bution,  /(d)  At;  particles  in  element  Ad.  As  particles  penetrate  the  plate,  the  simulation 
provides  the  velocity  of  each  particle  passing  the  last  grid  point  per  unit  time,  i.e.  the 
flux  distribution,  d/(d)Ad  particles  in  the  element  Ad  per  unit  time  At.  For  each  passing 
particle,  indexed  with  i,  then  d,/,(d)Ad  =  1.  If  we  let  Vi=vifl(v)Av,  then  the  particle  flux, 
in  numbers  per  At,  hitting  the  plate  is 

Fo  =  ]Tv;.  (14) 

i 

Similarly  the  effective  density  at  the  plate  is 

no  =  Y.  K/uj.  (15) 

« 

Hence  the  mean  speed  of  these  particles  is 

(v)  =  Fo/no.  (16) 

We  calculate  the  average  kinetic  energy  that  a  particle  has  at  the  plate  with 

*  -££«*!•  <I7> 

Finally  the  kinetic  energy  flux  is  determined  with 

Qo  =  yI>.2^  (18) 

t 

To  reduce  noise  in  these  diagnostics,  we  evaluate  the  spatially-averaged  plasma  period  and 
then  average  the  above  values  over  that  plasma  period. 

TV.  Simulation  Set-up 

We  use  the  same  model  as  that  shown  in  the  previous  report  (QPR  III, IV  1984). 


Specifically  the  input  used  follows  in  Table  1. 


Table  1.  Simulation  Set-Up 


System  Parameters 

Injection  Parameters 

Initially  Empty 

Ion  Flux=500. 

Open  External  Circuit 

Electron  Flux=500. 

Length=2. 

Electron  Thermal  Velocity =1 

Grid  Number=128 

Ion  Thermal  Velocity=0.158 

At=3.906xl0“3 

Electron  Temperature=0.05 

Particle  Charge=0.05 

Electron  Mass=0.05 

Ion  Mass=2. 

Ion  Temperature=0.05 

V.  Comparison  of  Theory  and  Simulation 

To  determine  the  theoretical  values  for  kinetic  energy  flux,  we  inserted  the  potential 
calculated  in  the  simulation  into  Eqs.  (10)-(12).  The  magnitude  of  the  plate  potential 
measured  from  the  simulation  exceeds  that  predicted  in  Eq.  (13)  because  the  injected  ions 
enter  the  system  with  Tt  =  T  rather  than  T,  T.  The  first  results  for  each  variable  listed 
in  Table  2  is  that  which  theory  predicts.  The  second  result  is  output  from  PDW1.  The 
range  indicated  for  simulation  results  show  the  amplitude  of  oscillation.  All  values  from 
simulation  are  measured  after  the  total  number  of  particles  in  the  simulation  is  unchanging. 

We  present  in  Fig.  2  an  example  of  the  history  of  kinetic  energy  flux  to  the  plate 
for  both  ion  streams.  Here  R  =  40%  for  a  flux  history  that  spans  the  last  fourth  of  the 
simulation  in  time.  The  small  dip  at  time  =  87  occurs  because  the  ion-ion  two-stream 
instability  causes  the  ion  streams  to  interact  and  whirl.  (See  the  previous  report.)  This 
results  in  an  effective  reduction  in  kinetic  energy  flux  because  of  the  drop  in  density  from 
the  ion  hole  created  in  phase  space  which  hits  the  plate.  Figure  3  shows  the  ion  phase  space 
profile  at  six  times  during  the  simulation. 

VI.  Conclusions 


All  of  the  assumptions  used  in  the  derivation  of  the  time-independent  theory  of  ki¬ 
netic  energy  flux  to  a  boundary  with  ion  reflection  are  verified  by  our  simulation  using 


PDWl.  Even  with  the  onset  of  two-streaming,  the  average  value  of  the  transport  results 
is  unaffected. 

We  determine  that  the  electrons  do  indeed  carry  an  average  energy  of  IkT  to  the  plate 
as  verified  in  Table  2  by  the  results  for  Qe o.  However  with  the  diagnostic  of  average  particle 
energy  at  the  plate,  we  see  that  PDWl  consistently  generates  £eo~kT/2  rather  than  kT. 
In  addition  PDWl  generates  for  all  three  cases  an  electron  temperature  profile  that  drops 
from  kT  at  the  source  to  kT/ 2  at  the  plate.  Thus  the  history  of  electron  kinetic  energy 
at  the  plate  indicates  that  £eo  equals  the  electron  temperature  measured  at  the  plate.  We 
point  out  in  Table  2  that  Qi 0  is  also  accurately  predicted.  This  verifies  the  assumption 
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that  the  ion  energy  carried  from  the  bulk  plasma  to  the  boundary  is  kT  —  e<f> 0.  Again 


the  value  £,o  specifically  calculated  by  PDWl  is  overpredicted  by  theory.  If  we  use  the 
computational  results  for  £&  along  with  Eq.  (11),  then  we  find  £ioo  is  0.84  kT  for  J2  =  0% 
and  0.88  kT  for  R  =  20%  and  40%.  Consequently  we  ask  the  following  question  with  the 


hope  that  in  the  near  future  we  can  answer  it.  Is  the  average  energy  that  a  particle  carries 
from  the  bulk  plasma  to  the  boundary  that  same  as  the  average  kinetic  energy  that  a 


particle  has  when  hitting  the  plate? 
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Theory  and  Simulation  Results 


Variable 

Reference 
Eq.  () 

R  0% 

20% 

40% 

kT 

(13) 

0.93 

1.33 

1.77 

— 

1.04  ±  0.10 

1.40  ±  0.11 

1.82  ±  0.15 

n«o 

* 

627 

627 

627 

(15) 

630  ±  90 

630  ±  110 

630  ±  180 

Feo 

— 

500 

500 

500 

(14) 

500  ±  20 

500  ±  20 

490  ±  30 

Fro 

(9) 

0 

-125 

-333 

(14) 

0 

-125  ±  7 

-330  ±  25 

£e0 

(3) 

0.05 

0.05 

0.05 

(17) 

0.026  ±  0.004 

0.025  ±  0.005 

0.025  ±  0.005 

£iO 

** 

0.102 

0.120 

0.141 

(17) 

0.094  ±  0.007 

0.114  db  0.004 

0.135  ±  0.006 

Qe 0 

(4) 

25 

25 

25 

(18) 

25  ±  1 

25  ±  1 

25  ±  2 

QiO 

(11) 

51 

75 

117 

(18) 

51  ±  2 

75  ±  4 

117  ±  7 

QrO 

(12) 

0 

-15 

-47 

(18) 

0 

-15  ±  3 

-47  ±  5 

neo  =  Fe0(2fcT/7rm)  1^2  from  Eqs.  (2)  and  (5). 
**  Use  Eq.  (10)  with  £tao  =kT. 


Table  2.  Comparison  of  theory  and  simulation  using 

the  values  presented  in  Table  1.  The  first  value  for 

each  variable  is  that  from  theory;  the  second  is  from  PDWl. 


Figure  3.  Profiles  in  ion  phase  space  at  various  times  for  R  =  40%.  Source  at  x  =  0 
Plate  at  x  =  2.  Ion  thermal  velocity  at  x  =  0  is  0.158.  Primary  ions  are  denoted  by 
reflected  ions  by  “+” . 
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B.  Planar  Magnetron  Discharges 

Amy  Wendt  (Prof.  M.  A-  Lieberman,  Dr.  Herman  Meuth) 

Perry  Gray  (ftof.  C.  K.  Birdsall) 

(1)  Progress  on  Experimental  Apparatus 

Construction  and  testing  of  the  experimental  apparatus  has  been  completed,  and  it  is  now  in 
use.  No  major  design  changes  from  that  already  described  in  previous  reports  were  necessary.  Elec¬ 
trodes  have  been  mounted  an  pistons  which  cany  two  teflon  o-rings  that  seal  the  piston  radially 
aggtn*  a  cylindrical  outer  chamber.  A  diffusion  pump  evacuates  the  main  vacuum  region,  and  a 
mechanical  p™p  differentially  pumps  between  the  two  s^ak  an  earh  piston.  "Ibis  prevents  air  from 
lairing  in»fi  the  mam  vacuum  regjon  when  the  pistons  are  moved  diming  operation.  This  design  has 
proven  successful;  the  system  has  a  base  pressure  of  less  than  10~s  ton. 

The  electromagnet  for  planar  magnetron  operation  of  the  experiment  has.  been  designed  and 
built.  The  design  was  completed  with  the  use  of  POISSON,  a  two  dimensional  magnet  meddling 
code  available  to  us  through  National  Magnetic  Fusion  Energy  Computing  Center.  Given  the 
geometry  of  a  magnet  design,  the  code  maps  out  magnetic  held  fines  and  contours  of  constant  mag¬ 
netic  field.  This  allowed  us  to  choose  a  configuration  of  magnetic  windings  and  iron  core  which 
would  fit  the  space  available  and  which  would  optimize  the  narrating  of  electrons  in  the  magnetic 
field,  hi  order  to  retain  same  flexibility  in  the  geometry  of  the  magnetic  field,  removable  pole  pieces 
were  used  an  the  iron  core  which  may  be  replaced  by  pole  pieces  of  different  shapes.  The  magnet 
has  been  wound  and  tests  have  been  made  of  the  water  cooling  and  dectrical  systems.  The  magnet 
is  now  set  up  an  a  workbench  so  that  its  magnetic  field  can  be  measured.  A  poor  man’s  optical 
bench  has  been  devised  so  that  a  gauss-meter  can  be  used  to  make  spatially  resolved  measurements 
of  radial  and  azimuthal  fields.  B,  was  measured  as  a  function  of  r  far  various  values  of  z  with 
results  very  similar  to  Rg.  3  in  QPR  3  and  4  1984,  page  S3;  B,  »1000  gauss  at  the  cathode  plate  was 
obtained  far  magnet  current  of  2900  A.  The  magnet  will  be  installed  as  soon  as  testing  is  completed. 

At  the  same  tune,  work  has  been  done  on  characterizing  the  dc  glow  discharge.  The  automatic 
pressure  control  system  has  been  set  up  and  works  wdL  Several  Iangmuir  probes  have  been  made 
and  are  mounted  an  the  diagnostic  parts.  The  probes  have  sliding  seals  so  the  radial  position  of  the 
probes  can  be  varied.  A  single  Langmuir  probe  has  been  made  out  of  0.005  in.  diameter  platinum 
wire  encased  in  alumina  tubing  with  a  0.014  in.  sphere  exposed  to  the  plasma.  The  spherical  probe 
tip  was  made  by  heating  the  end  of  the  platinum  wire  with  a  torch.  tW>  double  probes  have  also 
been  made.  The  two  probes  which  make  up  the  double  probe  are  biased  with  respect  to  each  other 
but  float  with  respect  to  the  chamber.  One  of  the  double  probes  was  made  with  platinum  wires 
exactly  as  described  above  far  the  single  probe.  The  other  double  probe  was  made  using  a  pair  of 
paralleled  platimim  disks  0.313  in.  in  diameter.  The  back  sides  of  the  discs  were  coated  with  epoxy 
resin  so  that  they  are  nonconductive.  The  increased  area  of  this  double  probe  allows  for  an  increased 
probe  current  so  that  measurements  can  be  made  even  at  low  density.  Circuitry  is  available  to  evalu¬ 
ate  the  probe  current-voltage  characteristic,  using  either  a  dc  bias  far  point-to-point  measurement  or 
a  swept  bias  far  continuous  measurement.  The  probe  has  been  used  to  measure  electron  density  in 
the  gjow  discharge  and  will  also  be  used  far  measuring  the  electron  temperature. 

We  have  also  studied  the  current-voltage  characteristic  of  the  discharge  at  various  pressures. 
We  have  observed  that  the  current  drawn  from  our  constant  voltage  supply  has  a  tendency  to  be 
unsteady  both  on  a  short  time  scale  (corresponding  to  flashes  of  light  from  the  discharge)  and  on  a 
longer  time  scale  (slow  drifts  of  both  increasing  arid  decreasing  current).  To  remedy  this  problem, 
we  are  increasing  the  voltage  dropping  resistors  which  are  in  series  with  the  glow  discharge  so  that 
effectively  the  discharge  will  see  a  constant  current  supply. 

(2)  Modelling  and  Simulation 

The  amulatian  model  we  are  working  with  far  the  planar  magnetron  has  been  described  in 
figure  6  page  56  in  the  previous  QPR  (3  and  4,  1984).  We  are  currently  implementing  modelling  of 


ionization  by  volume  creation  cl  electron-ion  pairs  with  provisions  to  make  this  a  function  of  local 
variables  (e.g.,  local  density),  allowing  the  magnetic  field  to  vary  with  distance  from  the  cathode, 
and  inserting  additional  transport  diagnostics. 

Electron-neutral  collisions  are  modelled  with  electrons  scattering  off  a  fixed  background  of 
neutrals  as  if  hard  and  massive  spheres  (elastic  coflisians).  At  a  given  time  step,  the  number  of  cofli- 
aons  that  occur  is  detrnrnnrd  through  an  inverse  cumulative  Poisson  distribution  so  as  to  be  random 
in  time,  about  a  given  mean  value.  Then  particles  to  be  scattered  are  chosen  randomly  from  the  par- 
tide  population;  a  particle  may  be  scattered  more  than  once  in  a  given  time  step.  Scattering  angles 
are  chosen  using  the  cosine  probability  far  hard  sphere  collisions  and  the  phase  angle  is  chosen  uni¬ 
formly  randomly  in  the  interval  0  to  2n.  The  computer  routines  to  model  coflisions  in  this  manner 
have  been  written  and  are  being  tested  to  see  whether  they  produce  the  expected  statistics  far  the 
angular  spread  in  time  of  an  initially  monrwenergeiir  distribution  of  halligtir  particles  (a  beam).  The 
free  parameter  for  this  modd  is  the  density  of  the  background  neutrals.  This  model  does  not  prop¬ 
erty  handle  energy  exchange  between  species  or  ionization  and  excitation  of  neutrals.  However,  we 
believe  it  is  a  good  first  step  far  the  effect  of  collisions  an  the  plasma  behavior  in  the  regime  where 
electron-neutral  coflisians  exceed  Coulomb  collisions. 


C.  Partide  Simulations  of  the  Low  a  Pierce  Diode 

T.  L.  Crystal,  S.  Kuhn 

This  work  was  completed  with  publication  of  an  aitide  with  the  above  title  and  authors  in 
Physics  of  Fluid,  28,  pp.  2116-2124,  July  1985. 


D.  Thermionic  Emission  I-V  Knee  Paradox 
Perry  Gray,  Win.  S.  Lawson  (Prof.  C.  K.  Hirdsafl) 

(1)  Introduction.  Electron  emissian  current  density  in  a  planar  diode  with  a  hot  dispenser 
cathode  (thermionic  emission)  is  supposedly  wdl  known.  However,  R.  T.  Longo1  proposed  that  the 
net  current  density  from  a  hot  cathode  was  a  particular  mixture  of  space-charge  and  temperature 
(united  values.  He  proposed  the  empirical  fit 


’  limited 


(i) 


where  the  space-charge  limited  current  density  is  the  usual  Child’s  Law  (1911) 


J5C  =KxV3/i/d2 


(2) 


and  the  temperature-limited  current  density  is  the  usual  Schottky  Law  (1914) 


Jtl 


=  J  wtmrted 


^V7/rl 


P) 


where  V  is  the  anode  voltage  and  d  is  the  cathode-anode  spacing.  Longer2  showed  remarkable  agree¬ 
ment,  reproduced  here  as  Figure  1,  with  Eq.  (1),  with  no  sharp  knee  juncture  of  Eqs.  (2,  3);  he 
labels  this  as  “Langmuir  Theory”  (perhaps  a  trifle  unfair). 
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SECTION  EH:  CODE  DEVELOPMENT  AND  SOFTWARE  DISTRIBUTION 


Codes  winch  we  have  developed  or  modified  here  are,  generally  speaking,  available  far  the  ask¬ 
ing.  We  are  making  arrangements  with  our  Industrial  Liaison  Program  to  handle  distribution,  at 
cost,  of  listings,  tapes  anri/nr  diskettes.  In  a  later  report  we  will  provide  a  list  of  codes  and  related 
material  available  and  how  to  obtain  such.  If  you  use  those,  we  simply  ask  that  you  acknowledge  the 
source  and  also  send  us  your  results  as  published  in  reports  or  journals. 

(Not  all  of  the  codes  mentioned  in  our  reports  are  available  far  distribution;  same  are  single¬ 
purpose,  single- user,  too  stylized  far  outside  use.) 

A.  POLY:  A  Hybrid  Scheme  for  the  Solution  of  the  Vlasov  Equation 

K.  Theilhabcr 


I.  INTRODUCTION. 

In  this  report,  we  present  a  hybrid  numerical  scheme  for  the  simulation  of 
plasmas  where  one  component  can  be  modelled  by  linearised  fluid  equations,  the 
remainder  of  the  plasma  being  described  by  the  full  Vlasov  equation.  For  both 
plasma  components,  numerical  solutions  are  obtained  through  finite-difference 
methods. 

For  instance,  the  hybrid  scheme  is  well  adapted  to  the  study  of  the  initial 
phases  of  the  weak  warm  beam  instability,  where  a  weak  beam  interacts  with 
a  massive  bulk  plasma.  In  this  situation,  it  is  essential  to  retain  all  of  the 
phase-space  dynamics  of  the  beam.  On  the  other  hand,  the  bulk  plasma  merely 
sustains  plasma  oscillations,  and  it  is  appropriate  to  model  it  as  a  linear  fluid, 
thereby  saving  a  very  considerable  amount  of  mesh  size  and  computational  time. 

In  what  follows,  we  shall  present  normalized  equations  for  the  beam-plasma 
problem,  and  then  sketch  our  methods  of  solution.  It  should  be  kept  in  mind 
that  the  formulation  of  the  problem  can  be  generalized  in  a  straighforward 
fashion:  thus,  the  "beam"  can  just  as  well  represent  a  warm  but  stationary 
plasma,  and  its  density  need  not  be  small  compared  to  the  density  of  the  bulk 
plasma;  the  model  for  the  bulk  plasma  can  be  modified,  so  as  to  include  thermal 
dispersion  or  damping,  effects  which  we  are  ignoring  in  the  present  formulation 
of  the  scheme. 

II.  THE  PHYSICAL  PROBLEM  AND  ITS  NORMALIZATIONS. 

Our  plasma  is  one-dimensional  and  consists  of  one  species,  electrons  moving 
against  an  immobile  neutralizing  background  of  ions.  A  qualitative  sketch  of 
the  initial,  unperturbed  distribution  function  is  shown  in  Fig.(l),  in  which  the 
curve  representing  the  bulk  should  be  thought  of  as  a  very  narrow,  cold  beam 
at  v  =  0.  Because  our  code  is  aimed  at  exploring  the  initial  turbulence  created 
in  the  beam-plasma  instability,  we  can  think  of  the  beam  distribution  function 
as  a  tunable  source  of  free  energy. 

Let  the  total,  averaged  spatial  density  of  the  electrons  be  no,  with  this 
density  split  between  the  bulk,  no p,  and  the  beam,  nos,  nop  +  nog  =  no-  In 
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the  weak  warm  beam  instability,  the  beam  is  a  small  perturbation  to  the  bulk 
so  that  nos  ^  no p,  but  this  condition  is  not  assumed  in  the  numerical  scheme. 

We  model  the  bulk  evolution  by  the  coupled  linearised  continuity  and  mo¬ 
mentum  equations: 


-nP  =  -nos—  us,  (l) 

S4- <2> 

where  rip  and  up  are  the  linearised  fluid  variables.  It  is  straightforward  to  add 
thermal  dispersion  or  damping  terms  to  these  equations,  but  for  simplicity  we 
shall  keep  the  bulk  plasma  as  cold  and  collisionless. 

The  beam  is  described  by  a  distribution  function,  fB  (*,  v,  t)  which  evolves 
according  to  the  Vlasov  equation: 

TtfB+vTxfB+mTvfB^Q'  M 

with  the  normalisation  of  the  unperturbed  distribution  function: 

f  dvfB(v)-noB,  (4) 

J  —  OO 

The  electric  field  E(x,  t)  is  determined  self-consistently  from  Poisson’s  equation: 

d  q  r 

—E(x,t)  =  —  t^p +hB(x,t)+  I  dvfB(x,v,t)  -  no  ,  (5) 

where  we  have  subtracted  from  the  electron-charge  source  terms  the  neutralizing 
charge  of  the  immobile  background  ions. 

We  proceed  to  normalise  eqs.(l-  5).  We  already  have  a  characteristic  time- 
scale,  determined  by  the  plasma  frequency: 


but,  because  the  bulk  is  cold,  we  cannot  choose  its  thermal  velocity  as  a  refer¬ 
ence  velocity,  nor  the  Debye  length  =  vr/uP  asa  reference  length.  Rather, 
we  introduce  a  reference  velocity  vB  which  will  remain  arbitrary  for  the  mo¬ 
ment,  but  which  eventually  will  be  tied  to  the  width  of  the  beam.  The  reference 
length  is  then  X  =  vB/u>p,  and  we  define  the  normalised  variables  as: 


t'  =  Wpt, 
x'  =  x/Xr, 
v'  =  v/vr. 


(7) 

(8) 
(9) 


O’ 

II 

(10) 

mujpVR 

u'p  =  up/vr, 

(ID 

n'p  =  n/»/nop, 

(12) 

/'  =  —  is- 

(13) 

tlOB 

Eqs.(l-  5)  then  become: 

d  ,  d  , 

avnp  =  ~^p' 

(14) 

T7“'p  =  &■ 
dt‘  p 

(15) 

a?7  at'7  av'7 

(16) 

1  dv’Mv’)  =  1. 

(17) 

=  if,  n',,  +  (I  /'  do-  l)  . 

(18) 

where  we  have  defined  the  important  parameters, 

p  _  *0fl 

H-B  » 

(19) 

^0 

Rp  =  m. 

no 

(20) 

which  gauge  the  relative  bulk  and  beam  densities,  the  total  density  being  nor* 
malised  to  1,  with  Rp  +  Rb  =  1.  For  convenience,  we  shall  drop  all  primes  in 

referring  to  the  normalised  equations. 

For  eqs.(14-  18),  we  can  define  the  normalized  energies  of  the  electric  field, 

of  the  bulk  plasma  and  of  the  beam  plasma  particles.  These  are, 

respectively: 
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and  with  these  definitions,  conservation  of  the  total  energy  of  the  system  follows: 


—  u»tot  =  +  «"b)  =  0, 


(24) 


m.  BOUNDARY  CONDITIONS. 

The  variables  np  and  up  are  defined  on  the  one  dimensional  spatial  mesh 
x i  =  (t  —  l)Ax,t  =  1,2,3 The  beam  distribution  function  f(x,v,t) 
requires  of  course  an  additional  mesh  in  v-space,  v}  —  +  (j  —  l)Av,j  = 

1, 2, 3, . . .  jVv,  where  ojrv  =  Umax,  and  /  *»  defined  on  the  JVX  X  Ny  two- 
dimensional  mesh.  The  spatial  boundary  conditions  are  chosen  to  be  periodic, 
while  in  v-space  we  require  the  distribution  function  to  be  sero  at  the  extreme 
boundaries,  v  =  Uj^  and  v  =  %ix-  This  latter  requirement  is  not  physical 
to  the  extent  that  the  distribution  function  should  really  extend  to  arbitrarily 
large  velocities.  In  operation,  the  Vlasov  solver  will  tend  to  transport  some 
density  beyond  the  maximum  velocity  bounds,  and  in  so  doing  will  "forget”  the 
profile  of  what  was  loet  to  higher  velocities,  resulting  in  a  density  (and  energy) 
loss  over  time.  However,  by  choosing  and  t>max  somewhat  larger  than 

the  beam  width  and  thus  insuring  that  /  is  very  small  at  the  boundaries,  these 
losses  can  be  made  extremely  small,  thus  insuring  consistency  of  the  assumption 
of  stricly  sero  density  at  the  edges.  The  exceptions  are  those  when  very  strong 
trapping  occurs,  in  which  case  the  trapping  widths  may  extend  to  beyond  the 
boundaries  and  density  lost  to  high  velocities. 


iv.  THE  POISSON  SOLVE E- 

We  solve  Poisson’s  equation: 


—E  =  n, 
ox 


in  two  steps,  through  a  potential  <j>'. 


a3 

a*2' 


'  =  -n, 


E  =  — —  <4. 

with  the  finite-difference  form  of  these  equations: 

4i+ 1  —  +4i- 1  =  —  A  x7iti 


Ei  =  ~2Ax^‘+l  ~  ^*-1^ 

With  a  Fast  Fourier  Transformation,  the  solution  for  E  is: 

.Ax  /ibAxN  ... 
E(k)  =  -t— cot  J  n(k), 

where  k  is  the  discrete  wave-vector. 


(25) 

(26) 
(27) 

(28) 

(29) 

(30) 


V.  PUSHING  THE  BULK  PLASMA. 

In  advancing  eqs.(14-  15),  we  use  a  finite-difference  scheme  of  the  ’’leapfrog* 
type,  with  the  density  defined  at  the  same  time  as  the  electric  field,  and  the 
velocity  at  ±At/2  time  intervals  before  and  after.  Dropping  the  P  subscripts 
on  the  variables,  we  have: 


At  ,  (n) 
•+ 1 

(™)  _ 


(«+l) 

u-  '  —  u 


(31) 

(32) 


where  the  superscripts  indicate  the  time  positioning.  We  solve  eqs.(31)  and 
(32)  in  Fourier  space,  with  the  FFT  equations: 


n(»+l/3)  (jfc)  _  n(n-l/3)  (jfc)  =  _t  £i  8in(fcA*)  «<">  (Jfe),  (33) 

u(»+l)  (jfc)  _  u<»)  (Jfc)  =  £<»+l/2)  (jfc).  (34) 

While  we  have  assumed  so  far  that  the  bulk  plasma  is  perfectly  cold,  the 
inclusion  of  a  thermal  dispersion  term  in  eqs.(33-  34)  should  be  straightforward, 
and  will  be  attempted  in  future  work.  The  overall  accuracy  of  eqs.(33-  34)  is 
0(  At2,  As3). 

VI.  PUSHING  THE  BEAM  PLASMA. 

The  beam  plasma  distribution  function  is  evolved  synchroneously  with  the 
bulk  equations  (31*  32),  using  a  split-step  method  (C.Z.  Cheng  and  G.  Knorr) 
which  involves  transport  in  only  one  dimension,  z  or  o,  at  a  time.  Had  we  a 
perfect  one-dimensional  transport  scheme,  the  code  would  proceed  as  follows: 


/(«+l/3)  (S)  „)  =  f(n)  (x  _  uAt)  „)  (35) 

/(n+1>(z,v)  =  /(«+i/2)(x,„  _  E{n+l/2)(x)At).  (36) 

where  i?<"+l/2>(x)  is  determined  self-conaistently  with  the  bulk  and  beam  densi¬ 
ties  njr+1^3*  and  /<n+1/2)(X)  v).  Such  a  procedure  can  be  shown  to  be  accurate 
to  order  0(At3),  the  main  difficulty  residing  in  finding  a  comparably  accurate 
transport  scheme.  Such  a  scheme  is  the  'flux-corrected  transport*  (FCT)  of 
Boris  and  Book(1),  which  incorporates  exceptionally  low  diffusivity  over  a  large 
number  of  time  steps. 

We  shall  not  describe  the  FCT  scheme  in  detail,  but  note  the  salient  fea¬ 
tures  of  the  version  used  here.  At  the  heart  of  FCT  is  a  'smart”,  nonlinear 
filter  which  is  interposed  between  a  diffusive  and  an  antidiffusive  step.  The 
nonlinear  filter  recognizes  which  features  of  the  transported  function  have  to  be 
diffused  for  numerical  stability,  but  preserves  to  the  maximum  extent  physically 
meaningful  features  such  as  shock  fronts.  Because  of  this,  FCT  is  well  adapted 
to  the  'mixing'  character  of  the  Vlasov  equation,  where  regions  of  different 
densities  may  coexist  on  short  space  and  velocity  scales. 

In  the  terminology  of  (*)  we  are  using  a  'low  phase  error',  explicit,  'phoeni- 
cal”  FCT  scheme,  which  solves  the  one-dimensional  transport  equation: 

T>r  +  Vk“  =  °'  (JT| 

In  the  notation  of  (1),  the  diffusion  and  antidiffusion  coefficients  are,  respec¬ 
tively: 


•  V  -*■  V 

VaV-V-V 


1  lj 

(38) 

1  1 2 
^=6'6e- 

(39) 

where  e  =  VAt/A£  measures  the  displacement  at  each  time  step  relative  to  a 
mesh  sise.  Transport  in  x  and  v  are  not  handled  in  exactly  the  same  way,  as 
in  the  first  case  periodic  and  in  the  latter  sero-value  boundary  conditions  are 
imposed.  In  particular,  as  discussed  in  Section  III,  some  particle  loss  to  the 
edges  in  velocity  space  is  inevitable  but  can  be  made  very  small. 

VTL  ENERGY  CONSERVATION. 

In  the  absence  of  the  beam  plasma,  the  bulk  equations  (31-  32),  coupled 
to  Poisson’s  equation  (30),  have  an  exactly  conserved  quantity  which  reduces 
to  the  sum  of  kinetic  and  electric-field  energies  in  the  limit  Ax,  At  — *  0.  If  we 
define  all  quantities  at  the  (n  +  1/2)  time-step,  then  the  electric-field  energy  is: 

and  the  corresponding  bulk  kinetic  energy  is: 


and  with  Rg  /  0,  the  quantity: 


+  u,<r+l/2)  +  ^n+1/2)  (44) 

is  conserved  to  £?(At2)  at  each  time  step.  These  quantities  correspond  to  the 
continuum  quantities  of  eqs.(21-  23),  and  eq.(44)  to  the  exact,  physical  energy 
conservation  relation  (24). 

VIU.  LANDAU  DAMPING  ON  A  MAXWELLIAN. 

For  the  initial  testing  of  the  code,  and  to  motivate  the  rational  development 
of  its  diagnostics,  we  considered  a  series  of  small-amplitude  numerical  experi¬ 
ments  in  which  the  accuracy  of  the  integration  scheme  could  be  confronted,  at 
least  in  the  domain  of  linear  theory. 

We  started  with  a  series  of  experiments  on  linear  Landau  damping,  for 
which  we  ignored  the  bulk  package  by  setting  Rp  —  0,  and  for  which  the 
’beam”  was  no  longer  a  beam,  but  a  Maxwellian  centered  at  v  =  0.  For  these 
experiments,  the  initial,  unperturbed  distribution  function  is  given  by: 

(45) 

and  we  considered  the  decay  of  a  small,  initial  perturbation: 

/(x,  «,t  =  0)  =  /o(o)  (1  +  6  cos kbx)  (46) 

where  6  <  1.  Note  that  in  our  normalised  system,  where  uiP  =  l,  the  Debye 
length  for  this  thermal  plasma  is  simply  At%.  Starting  from  (46),  an  exact 
analytic  solution  of  the  linearised  Vlasov- Poisson  system  yields: 

S(z,  0  =  Asin(fct2)coa(wRi)e-,«  +  Re£e“‘*  /  ^[ZZ)  dv<  (4?) 

where  w  =  ojr  +  t'7  is  the  exact  solution  of  the  plasma  dispersion  relation: 

((k,u,)  =  l  +  y  -^A-dv  =  0,  (48) 

which,  for  the  Maxwellian  distribution  function  becomes: 

1  +  *3 Alls2  *  1  +  )=0* 

where  f  =  w/\/2kAvt.  In  deriving  eq.(47),  we  did  assume  that  w/j/fct,  >■  Av&. 
The  first  term  in  (47)  is  the  linear  fluid  response  of  the  plasma,  with  in  addition  a 


Landau  damping  term  exp(-'jt).  The  second  term  is  the  ballistic  contribution, 
which  for  the  Maxwellian  /0(t>)  decays  as: 


(50) 

where  the  phase-mixing  time  tpu  is  given  by: 

For  instance,  we  analysed  the  time-evolution  of  the  electric  field  energy 
with  the  initial  amplitude  b  =  0.001.  For  t  >  tpu  we  should  have: 

*»£(*)  =  tu£(0)e~2l'cos3(u>*t),  (52) 

and  thus  an  analysis  of  the  curve  of  wr  (t)  can  yield  the  numerical  values  of 
WR  and  7.  An  artisanal  but  reliable  technique  was  used:  the  real  part  of  the 
frequency  was  found  by  measuring,  by  hand,  the  spacing  of  the  seros  of  the 
field  energy,  while  the  damping  was  found  from  the  decay  of  the  envelope  of 
the  curve,  also  measured  by  hand.  In  fig.(2),  the  results  of  several  runs  are 
compared  to  the  exact  solution  of  eq.(49),  obtained  through  Newton’s  method. 
The  agreement  is  quite  good.  As  an  aside,  we  should  note  that  the  traditional 
electron-plasma  dispersion  relations: 


(51) 


ui  =  1  + 

Is 


(jJ 

y/2  ki,  A  vi,  ’ 


is  of  limited  validity,  k^Avi  <  0.1. 


(53) 

(54) 


DC.  CONCLUSION. 

The  predictions  of  the  hybrid  code  POLY  have  been  checked  in  the  case  of 
Landau  damping,  a  nontrivial  linear  phenomenon.  Furthermore,  in  connection 
with  the  study  of  beam-plasma  turbulence  (the  "turbulent  trapping  model*  ^ 

we  have  been  able  to  verify  that  the  code  yields 
results  for  small-scale,  nonlinear  plasma  behavior  which  are  consistent  with  the 
analytic  models.  We  are  thus  gaining  confidence  that  POLY  is  a  robust  numer¬ 
ical  tool,  which  should  have  a  wide  applicability  in  situations  where  detailed 
dynamics  are  required  for  only  a  limited  region  of  phase  space. 
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Numerical  Simulations  of  "Turbulent  Trapping"  and  the  Bump-on- Tail 
Ins tability .  K.  Theilhaber,  Electronics  Research  Laboratory,  University 
of  California,  Berkeley,  CA  94720;  D.  Pesrae  and  G.  Laval,  Centre  de 
Physique  Theorique,  Ecole  Poly technique,  91128  Palaiseau,  France. 


Theoretical  studies  (1)  strongly  suggest  that  in  the  presence  of 
sufficient  turbulence,  the  growth  rate  of  the  weak  warm-beam  instability 
is  significantly  modified  from  its  quasilinear  value,  and  that  a  con¬ 
current  modification  of  velocityrspace  diffusion  should  occur.  The 
criterion  is  approximately  (k  D)^  >  ,  where  D  is  the  quasilinear 

diffusion  coefficient,  and  the  quasilinear  growth  rate.  This  is  in 
fact  a  regime  always  attained  after  a  sufficient  number  of  e-foldings. 

We  numerically  study  the  old  paradigm  of  quasilinear  theory  by  starting 
the  instability  In  the  presence  of  an  initial  turbulence,  for  a  one-species, 
one-dimensional  configuration.  We  shall  present  results  for  growth- ra tes , 
diffusion,  and  the  correlation  functions  of  the  distribution  function. 

In  particular,  we  shall  show  the  links  to  the  "clump"  turbulence  of 
T.H.  Dupree  and  co-workers  (2). 
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(2)  T.H.  Dupree,  Phys.  of  Fluids,  15,  334  (1972). 


Simulations  of  the  Effect  of  ICRF  Waves  on  the 
Interchange  Instability* 


Niels  F.  Otani 

Electronics  Research  Laboratory 
University  of  California 
Berkeley,  CA  94720 


Recent  results  from  the  Phaedrus  experiment  at  Wisconsin  suggest  that  stabiliza¬ 
tion  of  the  interchange  mode  in  mirror  plasma  may  be  possible  using  RF  waves 
excited  externally  in  the  vicinity  of  the  ion-cyclotron  frequency.1  Two  simulation  codes 
have  been  developed  for  the  purpose  of  studying  the  effect  of  ICRF  waves  on  the 
linear  growth  rate  of  the  interchange  instability.  Both  codes  model  essentially  the 
same  system.  The  simulation  plane  is  oriented  perpendicular  to  the  background  mag¬ 
netic  field.  An  external  gravitational  field  is  fixed  antiparallel  to  the  density  gradients 
initially  present  in  the  system.  Full  dynamics  of  the  ions  are  followed  while  the  elec¬ 
trons  are  represented  as  a  cold  ExB  fluid.  The  plasma  is  assumed  quasineutral  and 
the  Darwin  approximation  is  made.  In  the  original  code,  no  significant  stabilization  is 
observed  for  long  interchange  wavelengths  (kL„^.\)  and  large  RF  fields 
(e2V|£|2/4m,2  «2S  =  0(1)).  The  parameters  required  for  tests  of  the  effects  on  short 
wavelength  modes  (/c£„»l)  excite  numerical  instabilities  which  appear  related  to  an 
insufficiently  accurate  representation  of  the  Lorentz  force.  A  new  code  employing  a 
different  algorithm  is  being  developed  to  bypass  this  problem.  Supporting  theory  will 
also  be  presented. 

‘Work  supported  by  DOE  Contract  No.  DE-AT03-76ET53064.  Computations  facilities 
provided  by  the  National  Magnetic  Fusion  Energy  Computer  Center. 

'•I.  R.  Ferron,  N.  Hershkowitz,  R.  A.  Breun,  S.  .  Golovato,  and  R.  Goulding  Phvs. 
Rev.  Lett.  51,  1955  0983). 


A  MODEL  OP  THE  PLASMA-SHEATH  REGION  INCLUDING, 
SECONDARY  ELECTRON  EMISSION  AND  ION  REFLECTION 


L.A.Schwager  and  C.K.Birdsall 

Plasma  Theory  and  Simulation  Group 
Electronics  Research  Laboratory 
University  of  California,  Berkeley  CA  94720 


A  low 'temperature  plasma  interacts  with  a  collector  plate 
primarily  through  the  emission  of  secondary  electrons  and 
the  reflection  of  plasma  ions.  These  phenomena  affect  the 
electrostatic  potential  drop  across  the  sheath  region  and 
hence  plasma  transport  to  the  plate.  The  computer  code  PDW1 
models  this  behavior.  The  code  is  a  one-dimensional,  time- 
dependent  electrostatic  particle  simulation.  Particle  ions 
and  electrons,  with  a  reduced  mass  ratio,  are  injected  con¬ 
tinuously  from  a  plasma  region  into  a  plate  where  ions  are 
reflected  or  cold  electrons  are  emitted.  A  parametric  study 
of  the  dependence  of  potential  drop  on  interaction  coeffi¬ 
cient  and  mass  ratio  is  presented.  We  have  also  developed 
theory  for  the  effect  on  potential  of  ion  reflection,  for 
the  cases  of  Ti=0  and  Ti=Te,  and  of  secondary  electron  emis¬ 
sion  for  Ti=Te.  The  simulation  results  for  ion  reflection 
approach  those  of  our  theory  for  the  cold  ion  case.  The  sim¬ 
ulation  results  with  secondary  electrons  compare  well  with 
previous  theory1  for  cold  ions  and  demonstrate  charge  sat¬ 
uration  as  predicted. 

^.D. Hobbs  and  J. A. Wesson,  Plasma  Physics,  Vol.9,  86  (1967) 

Work  supported  by  U .S .Department  of  Energy  under 
contract  DE-AT03-76ET53064. 

Computation  facilities  provided  by  NMFECC. 
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